ST465– Making Maps in R

Today we’ll go over how to make maps in R and how to work with spatial data sets. Some of the key R packages that we’ll work with are tmap, mapview for map making and sf, spData, spDataLarge for working with spatial data.

## R packages 
setwd(getwd())
library(rmarkdown)
library(tmap) ## Very commonly used -- to make static and interactive plots
library(mapview)
library(leaflet)
library(shiny)
library(tidyverse)
library(spDataLarge) ## Data sets -- use commands in notes to download
library(spData) ## Data sets
library(sf) ## simple features 

tmap

Let’s work with the world and coffee data sets. For each, identify the type of object it is and its coordinate reference system (CRS).

#data("world")
world
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 11
##    iso_a2 name_l~1 conti~2 regio~3 subre~4 type  area_~5     pop lifeExp gdpPe~6
##  * <chr>  <chr>    <chr>   <chr>   <chr>   <chr>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji     Oceania Oceania Melane~ Sove~  1.93e4  8.86e5    70.0   8222.
##  2 TZ     Tanzania Africa  Africa  Easter~ Sove~  9.33e5  5.22e7    64.2   2402.
##  3 EH     Western~ Africa  Africa  Northe~ Inde~  9.63e4 NA         NA       NA 
##  4 CA     Canada   North ~ Americ~ Northe~ Sove~  1.00e7  3.55e7    82.0  43079.
##  5 US     United ~ North ~ Americ~ Northe~ Coun~  9.51e6  3.19e8    78.8  51922.
##  6 KZ     Kazakhs~ Asia    Asia    Centra~ Sove~  2.73e6  1.73e7    71.6  23587.
##  7 UZ     Uzbekis~ Asia    Asia    Centra~ Sove~  4.61e5  3.08e7    71.0   5371.
##  8 PG     Papua N~ Oceania Oceania Melane~ Sove~  4.65e5  7.76e6    65.2   3709.
##  9 ID     Indones~ Asia    Asia    South-~ Sove~  1.82e6  2.55e8    68.9  10003.
## 10 AR     Argenti~ South ~ Americ~ South ~ Sove~  2.78e6  4.30e7    76.3  18798.
## # ... with 167 more rows, 1 more variable: geom <MULTIPOLYGON [arc_degree]>,
## #   and abbreviated variable names 1: name_long, 2: continent, 3: region_un,
## #   4: subregion, 5: area_km2, 6: gdpPercap
coffee_data
## # A tibble: 47 x 3
##    name_long                  coffee_production_2016 coffee_production_2017
##    <chr>                                       <int>                  <int>
##  1 "Angola"                                       NA                     NA
##  2 "Bolivia"                                       3                      4
##  3 "Brazil"                                     3277                   2786
##  4 "Burundi"                                      37                     38
##  5 "Cameroon"                                      8                      6
##  6 "Central African Republic"                     NA                     NA
##  7 "Congo, Dem. Rep. of"                           4                     12
##  8 "Colombia"                                   1330                   1169
##  9 "Costa Rica"                                   28                     32
## 10 "C\u00f4te d'Ivoire"                          114                    130
## # ... with 37 more rows
world_coffee = dplyr::left_join(world, coffee_data, by = "name_long")
world_coffee
## Simple feature collection with 177 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## # A tibble: 177 x 13
##    iso_a2 name_l~1 conti~2 regio~3 subre~4 type  area_~5     pop lifeExp gdpPe~6
##    <chr>  <chr>    <chr>   <chr>   <chr>   <chr>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 FJ     Fiji     Oceania Oceania Melane~ Sove~  1.93e4  8.86e5    70.0   8222.
##  2 TZ     Tanzania Africa  Africa  Easter~ Sove~  9.33e5  5.22e7    64.2   2402.
##  3 EH     Western~ Africa  Africa  Northe~ Inde~  9.63e4 NA         NA       NA 
##  4 CA     Canada   North ~ Americ~ Northe~ Sove~  1.00e7  3.55e7    82.0  43079.
##  5 US     United ~ North ~ Americ~ Northe~ Coun~  9.51e6  3.19e8    78.8  51922.
##  6 KZ     Kazakhs~ Asia    Asia    Centra~ Sove~  2.73e6  1.73e7    71.6  23587.
##  7 UZ     Uzbekis~ Asia    Asia    Centra~ Sove~  4.61e5  3.08e7    71.0   5371.
##  8 PG     Papua N~ Oceania Oceania Melane~ Sove~  4.65e5  7.76e6    65.2   3709.
##  9 ID     Indones~ Asia    Asia    South-~ Sove~  1.82e6  2.55e8    68.9  10003.
## 10 AR     Argenti~ South ~ Americ~ South ~ Sove~  2.78e6  4.30e7    76.3  18798.
## # ... with 167 more rows, 3 more variables: geom <MULTIPOLYGON [arc_degree]>,
## #   coffee_production_2016 <int>, coffee_production_2017 <int>, and abbreviated
## #   variable names 1: name_long, 2: continent, 3: region_un, 4: subregion,
## #   5: area_km2, 6: gdpPercap

We’ve created a new data set called world_coffee, what type of object is it and what is its CRS?

Let’s plot the spatial object using the R package tmap.

## tmap_mode "view" set to interactive viewing
tmap_mode("plot")
facets = c("coffee_production_2016", "coffee_production_2017")

tm_shape(world_coffee) + 
  tm_polygons(facets) + 
  tm_facets(nrow = 1, sync = TRUE)

Now let’s use the World data set (not that world and World are different data sets.)

data("World")
World 
## Simple feature collection with 177 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS:  WGS 84
## First 10 features:
##    iso_a3                   name           sovereignt               continent
## 1     AFG            Afghanistan          Afghanistan                    Asia
## 2     AGO                 Angola               Angola                  Africa
## 3     ALB                Albania              Albania                  Europe
## 4     ARE   United Arab Emirates United Arab Emirates                    Asia
## 5     ARG              Argentina            Argentina           South America
## 6     ARM                Armenia              Armenia                    Asia
## 7     ATA             Antarctica           Antarctica              Antarctica
## 8     ATF Fr. S. Antarctic Lands               France Seven seas (open ocean)
## 9     AUS              Australia            Australia                 Oceania
## 10    AUT                Austria              Austria                  Europe
##                   area  pop_est pop_est_dens                    economy
## 1    652860.000 [km^2] 28400000 4.350090e+01  7. Least developed region
## 2   1246700.000 [km^2] 12799293 1.026654e+01  7. Least developed region
## 3     27400.000 [km^2]  3639453 1.328268e+02       6. Developing region
## 4     71252.172 [km^2]  4798491 6.734519e+01       6. Developing region
## 5   2736690.000 [km^2] 40913584 1.495003e+01    5. Emerging region: G20
## 6     28470.000 [km^2]  2967004 1.042151e+02       6. Developing region
## 7  12259213.973 [km^2]     3802 3.101341e-04       6. Developing region
## 8      7257.455 [km^2]      140 1.929051e-02       6. Developing region
## 9   7682300.000 [km^2] 21262641 2.767744e+00 2. Developed region: nonG7
## 10    82523.000 [km^2]  8210281 9.949082e+01 2. Developed region: nonG7
##                 income_grp gdp_cap_est life_exp well_being footprint inequality
## 1            5. Low income    784.1549   59.668        3.8      0.79 0.42655744
## 2   3. Upper middle income   8617.6635       NA         NA        NA         NA
## 3   4. Lower middle income   5992.6588   77.347        5.5      2.21 0.16513372
## 4  2. High income: nonOECD  38407.9078       NA         NA        NA         NA
## 5   3. Upper middle income  14027.1261   75.927        6.5      3.14 0.16423830
## 6   4. Lower middle income   6326.2469   74.446        4.3      2.23 0.21664810
## 7  2. High income: nonOECD 200000.0000       NA         NA        NA         NA
## 8  2. High income: nonOECD 114285.7143       NA         NA        NA         NA
## 9     1. High income: OECD  37634.0832   82.052        7.2      9.31 0.08067825
## 10    1. High income: OECD  40132.6093   81.004        7.4      6.06 0.07129351
##         HPI                       geometry
## 1  20.22535 MULTIPOLYGON (((61.21082 35...
## 2        NA MULTIPOLYGON (((16.32653 -5...
## 3  36.76687 MULTIPOLYGON (((20.59025 41...
## 4        NA MULTIPOLYGON (((51.57952 24...
## 5  35.19024 MULTIPOLYGON (((-65.5 -55.2...
## 6  25.66642 MULTIPOLYGON (((43.58275 41...
## 7        NA MULTIPOLYGON (((-59.57209 -...
## 8        NA MULTIPOLYGON (((68.935 -48....
## 9  21.22897 MULTIPOLYGON (((145.398 -40...
## 10 30.47822 MULTIPOLYGON (((16.97967 48...

HPI stands for “Happy Planet Index” and provides a quantification of a country’s “happiness”. More information here: https://happyplanetindex.org Let’s check it out:

tmap_mode("view")

tm_shape(World) +
    tm_polygons("HPI")

Let’s look at the HPI along with economy type:

tmap_mode("view")
tm_shape(World) +
    tm_polygons(c("HPI", "economy")) +
    tm_facets(sync = TRUE, ncol = 2)

We continue to make maps, this time with the metro data set. What information is provided? What is the geometry type and what is the CRS?

## reading in data
data(metro)
metro 
## Simple feature collection with 436 features and 12 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -123.1193 ymin: -37.814 xmax: 174.7667 ymax: 60.16925
## Geodetic CRS:  WGS 84
## First 10 features:
##            name             name_long iso_a3 pop1950 pop1960 pop1970 pop1980
## 2         Kabul                 Kabul    AFG  170784  285352  471891  977824
## 8       Algiers El Djazair  (Algiers)    DZA  516450  871636 1281127 1621442
## 13       Luanda                Luanda    AGO  138413  219427  459225  771349
## 16 Buenos Aires          Buenos Aires    ARG 5097612 6597634 8104621 9422362
## 17      Cordoba               Cordoba    ARG  429249  605309  809794 1009521
## 25      Rosario               Rosario    ARG  554483  671349  816230  953491
## 32      Yerevan               Yerevan    ARM  341432  537759  778158 1041587
## 33     Adelaide              Adelaide    AUS  429277  571822  850168  971856
## 34     Brisbane              Brisbane    AUS  441718  602999  904777 1134833
## 37    Melbourne             Melbourne    AUS 1331966 1851220 2499109 2839019
##     pop1990  pop2000  pop2010  pop2020  pop2030                    geometry
## 2   1549320  2401109  3722320  5721697  8279607   POINT (69.17246 34.52889)
## 8   1797068  2140577  2432023  2835218  3404575     POINT (3.04197 36.7525)
## 13  1390240  2591388  4508434  6836849 10428756   POINT (13.23432 -8.83682)
## 16 10513284 12406780 14245871 15894307 16956491 POINT (-58.40037 -34.60508)
## 17  1200168  1347561  1459268  1562509  1718192  POINT (-64.18105 -31.4135)
## 25  1083819  1152387  1298073  1453814  1606993 POINT (-60.63932 -32.94682)
## 32  1174524  1111301  1065597  1023703  1057459     POINT (44.51462 40.182)
## 33  1081618  1141623  1217990  1320783  1505422  POINT (138.5986 -34.92866)
## 34  1381306  1666203  2033617  2388517  2721325  POINT (153.0281 -27.46794)
## 37  3154314  3460541  3951216  4500501  5070873    POINT (144.9633 -37.814)
tmap_mode("view")
tm_basemap("Stamen.Watercolor") +
tm_shape(metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")

You can also play around with the style in tmap.

tmap_style("classic")
## tmap style set to "classic"
## other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "watercolor"

tm_shape(World) +
    tm_polygons("HPI", legend.title = "Happy Planet Index")

We can also look through some data from New Zealand. As the others, what type of object is it and what is the CRS? How it different from others and why might that be the case?

## New Zealand
nz
## Simple feature collection with 16 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1090144 ymin: 4748537 xmax: 2089533 ymax: 6191874
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
## First 10 features:
##                 Name Island Land_area Population Median_income Sex_ratio
## 1          Northland  North 12500.561     175500         23400 0.9424532
## 2           Auckland  North  4941.573    1657200         29600 0.9442858
## 3            Waikato  North 23900.036     460100         27900 0.9520500
## 4      Bay of Plenty  North 12071.145     299900         26200 0.9280391
## 5           Gisborne  North  8385.827      48500         24400 0.9349734
## 6        Hawke's Bay  North 14137.524     164000         26100 0.9238375
## 7           Taranaki  North  7254.480     118000         29100 0.9569363
## 8  Manawatu-Wanganui  North 22220.608     234500         25000 0.9387734
## 9         Wellington  North  8048.553     513900         32700 0.9335524
## 10        West Coast  South 23245.456      32400         26900 1.0139072
##                              geom
## 1  MULTIPOLYGON (((1745493 600...
## 2  MULTIPOLYGON (((1803822 590...
## 3  MULTIPOLYGON (((1860345 585...
## 4  MULTIPOLYGON (((2049387 583...
## 5  MULTIPOLYGON (((2024489 567...
## 6  MULTIPOLYGON (((2024489 567...
## 7  MULTIPOLYGON (((1740438 571...
## 8  MULTIPOLYGON (((1866732 566...
## 9  MULTIPOLYGON (((1881590 548...
## 10 MULTIPOLYGON (((1557042 531...
legend_title = expression("Area (km"^2*")")
map_nza = tm_shape(nz) +
  tm_fill(col = "Land_area", title = legend_title) + tm_borders()
map_nza + tm_style("bw")
map_nza + tm_style("classic")
map_nza + tm_style("cobalt")
map_nza + tm_style("col_blind")

mapview

Let’s use mapview now, another common R package for constructing spatial maps. We’ll use the trails and franconia data set. What kind of objects are they and what is its CRS?

trails 
## Simple feature collection with 543 features and 3 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 528933.9 ymin: 5415754 xmax: 715816.7 ymax: 5597457
## Projected CRS: WGS 84 / UTM zone 32N
## First 10 features:
##                                                                   FGN    FKN
## 1                                                003756/Kunigundenweg 003756
## 2  004037/Jakobsweg (Almerswind-Coburg-Lichtenfels-Bamberg-Nuernberg) 004037
## 3                 005160/Steigerwaelder Jakobsweg (Bamberg-Uffenheim) 005160
## 4                                     022650/Sieben-Fluesse-Wanderweg 022650
## 5              005316/Rennsteigverein 1896 e.V. / Bamberger Rennsteig 005316
## 6                                      012029/Steigerwald-Panoramaweg 012029
## 7                            013633/Frankenwaldverein / Markgrafenweg 013633
## 8                                            023964/Rot-Main-Auen-Weg 023964
## 9                                                007443/Jean-Paul-Weg 007443
## 10                                               007443/Jean-Paul-Weg 007443
##       district                       geometry
## 1  Oberfranken MULTILINESTRING ((631670.9 ...
## 2  Oberfranken MULTILINESTRING ((636557.1 ...
## 3  Oberfranken MULTILINESTRING ((633443.8 ...
## 4  Oberfranken MULTILINESTRING ((631625.6 ...
## 5  Oberfranken MULTILINESTRING ((636894.6 ...
## 6  Oberfranken MULTILINESTRING ((634029.7 ...
## 7  Oberfranken MULTILINESTRING ((684491.3 ...
## 8  Oberfranken MULTILINESTRING ((686023 55...
## 9  Oberfranken MULTILINESTRING ((684861.6 ...
## 10 Oberfranken MULTILINESTRING ((684905.6 ...
franconia
## Simple feature collection with 37 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 8.975926 ymin: 48.8625 xmax: 12.27535 ymax: 50.56422
## Geodetic CRS:  WGS 84
## First 10 features:
##    NUTS_ID  SHAPE_AREA SHAPE_LEN CNTR_CODE                  NAME_ASCI
## 1    DE241 0.006736012 0.3926225        DE  Bamberg, Kreisfreie Stadt
## 2    DE242 0.008424469 0.6247263        DE Bayreuth, Kreisfreie Stadt
## 3    DE243 0.005982341 0.5185471        DE   Coburg, Kreisfreie Stadt
## 4    DE244 0.007329480 0.4569815        DE      Hof, Kreisfreie Stadt
## 5    DE245 0.146698316 3.4819699        DE         Bamberg, Landkreis
## 6    DE246 0.159489736 3.6242023        DE        Bayreuth, Landkreis
## 7    DE247 0.074698748 2.6954234        DE          Coburg, Landkreis
## 8    DE248 0.079746707 1.7712298        DE                  Forchheim
## 9    DE249 0.112934151 2.7544708        DE             Hof, Landkreis
## 10   DE24A 0.081960299 1.9393830        DE                    Kronach
##                          geometry    district
## 1  MULTIPOLYGON (((10.92581 49... Oberfranken
## 2  MULTIPOLYGON (((11.58157 49... Oberfranken
## 3  MULTIPOLYGON (((10.95355 50... Oberfranken
## 4  MULTIPOLYGON (((11.93067 50... Oberfranken
## 5  MULTIPOLYGON (((10.87615 50... Oberfranken
## 6  MULTIPOLYGON (((11.70657 50... Oberfranken
## 7  MULTIPOLYGON (((10.88654 50... Oberfranken
## 8  MULTIPOLYGON (((11.26376 49... Oberfranken
## 9  MULTIPOLYGON (((11.91988 50... Oberfranken
## 10 MULTIPOLYGON (((11.36979 50... Oberfranken

What does it look like if we only plot the trails data?

trails%>%
  mapview()

What about the franconia data set?

franconia%>%
  mapview()

We’d like to plot the trails onto the map, but they’re not the same CRS. In lecture, we saw that we can use a suite of functions starting with st_ to modify a spatial object. We’ll use a few here to change the crs and plot the data sets together:

## first, using st_transform() function we changed the CRS 
trails %>%
  st_transform(st_crs(franconia)) %>%
  st_intersection(franconia[franconia$district == "Oberfranken", ]) %>%
  st_collection_extract("LINE") %>%
  mapview(color = "red", lwd = 3, layer.name = "trails") +
  mapview(franconia, zcol = "district", burst = TRUE) 

Lastly, suppose you also wanted to know about any breweries along the trails. Let’s check out the breweries data set. What kind of object is it and what is its CRS?

breweries

What does a map of the breweries look like?

breweries%>%
  mapview()

Let’s combine them all:

trails %>%
  st_transform(st_crs(franconia)) %>%
  st_intersection(franconia[franconia$district == "Oberfranken", ]) %>%
  st_collection_extract("LINE") %>%
  mapview(color = "red", lwd = 3, layer.name = "trails") +
  mapview(franconia, zcol = "district", burst = TRUE) +
  breweries

tmap

What if we tried to do this in tmap?

tm_shape(franconia) + tm_fill("district") + 
tm_shape(trails) + tm_lines() + 
tm_shape(breweries) + tm_dots()